set.seed(12618)

2 c

Ejemplo ilustrativo de la verosimilitud de una Bernoulli vs nuestra Bernoulli truncada

2d

Graficar el sesgo del EMV en función del tamaño muestral n suponiendo que θ = 1/4.

tita <- 1/4
p_real <- 1/6 + 2/3 * tita
Nrep <- 10000
n_vals <- c(1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 60, 100, 200, 300, 400, 500, 1000, 1500, 2000)
sesgos <- numeric(length(n_vals))

esperanza_p <- function(n, p_real, k_esimo_momento = 1) {
  
  si_promedio_menor_o_igual_a_un_sexto <- function(n, p_real, k_esimo_momento = 1) {
    k <- floor(n / 6)
    return(((1/6)**k_esimo_momento) * pbinom(k, size = n, prob = p_real))
  }
  
  si_promedio_entre_un_sexto_y_cinco_sextos <- function(n, p_real, k_esimo_momento = 1) {
    lower <- floor(n / 6) + 1
    upper <- ceiling(5 * n / 6) - 1
    if (lower > upper) return(0)
    k_vals <- lower:upper
    probs <- dbinom(k_vals, size = n, prob = p_real)
    pesos <- (k_vals / n) ** k_esimo_momento
    
    return(sum(probs * pesos))
  }

  
  si_promedio_mayor_o_igual_a_cinco_sextos <- function(n, p_real, k_esimo_momento = 1) {
    lower <- ceiling(5 * n / 6)
    return(((5/6)**k_esimo_momento) * (1 - pbinom(lower - 1, size = n, prob = p_real)))
  }
  
  return(
    si_promedio_menor_o_igual_a_un_sexto(n, p_real, k_esimo_momento) +
    si_promedio_entre_un_sexto_y_cinco_sextos(n, p_real, k_esimo_momento) +
    si_promedio_mayor_o_igual_a_cinco_sextos(n, p_real, k_esimo_momento)
  )
}

esperanza_tita <- function(n, p_real) {
  return ((3/2) * (esperanza_p(n, p_real, 1) - 1/6))
}


for (i in seq_along(n_vals)) {
  n <- n_vals[i]
  estimadores <- numeric(Nrep)
  sesgos[i] <- esperanza_tita(n, p_real) - tita
}

Latex

Interactivo

El gráfico muestra el del estimador de máxima verosimilitud para distintos tamaños muestrales \(n\), con \(\theta = \frac{1}{4}\). Es más pronunciado cuando \(n\) es pequeño. A medida que \(n\) aumenta, el estimador tiende a concentrarse en torno a su valor esperado y el sesgo disminuye, tendiendo a cero.

Del gráfico podríamos deducir que el estimador resulta ser asintóticamente insesgado.

2 e

Graficar el ECM del EMV en función del tamaño muestral n suponiendo que θ = 1/4. Superponer un gráfico del ECM del EMV calculado utilizando el método clásico suponiendo que la proporción de estudiantes que se copió alguna vez en un examen y miente cuando se le hace la pregunta directa es π = 1/3.

proba_pi <- 1/3
p_y_clas <- (1 - proba_pi) * tita
Nrep <- 10000 
n_vals <- seq(2, 200, by = 2)
ecm_aleatorizado <- numeric(length(n_vals))
ecm_clas <- numeric(length(n_vals))


ecm_aleat <- function(n, p_real){
  sesgo_tita <- function(n, p_real){
    return (esperanza_tita(n, p_real) - tita) 
  }
  
  varianza_tita <- function(n, p_real){
    return (((3/2)**2) * (esperanza_p(n, p_real, k_esimo_momento = 2) - esperanza_p(n, p_real)**2))
  }
  
  return (varianza_tita(n, p_real) + sesgo_tita(n, p_real)**2)
}

for (i in seq_along(n_vals)) {
  n <- n_vals[i]
  
  ecm_aleatorizado[i] <- ecm_aleat(n, p_real)
  
  muestras_clas <- rbinom(Nrep, n, p_y_clas)
  promedios <- muestras_clas / n
  ecm_clas[i] <- mean((promedios - tita)^2)
}

Latex

Interactivo

El gráfico compara el del estimador de máxima verosimilitud (EMV) obtenido mediante el método de (curva azul), con el ECM del estimador (curva roja), suponiendo que la proporción verdadera de individuos con la característica de interés es \(\theta = \frac{1}{4}\) y que el porcentaje de estudiantes que mienten en el método clásico es \(\pi = \frac{1}{3}\).

2 f

Método numérico

p <- seq(1/6, 5/6, length.out = 20000)

l_diff <- -2 * 20 * log(p) - 2 * 80 * log(1 - p) + 
          2 * 20 * log(20 / 100) + 2 * 80 * log(80 / 100)
tita <- 3/2 * (p - 1/6)
corte <- qchisq(1-0.05, df = 1)

Latex

Interactivo

2g

Simular ambos experimentos en R suponiendo θ = 1/4 y p = 1/3, es decir, generar las respuestas de los encuestados por el método clásico y por el método de respuesta aleatorizada, para n = 10, n = 100 y n = 1000 y calcular el ECM empírico. Presentar los resultados en una tabla y describirlos, verificando que son compatibles con los resultados teóricos hallados en los ítems anteriores.

tita <- 1/4
proba_pi <- 1/3
p_y_aleat <- 1/6 + 2/3 * tita
p_y_clas <- (1 - proba_pi) * tita
n_vals <- c(10,100,1000,2000,3000)
Nrep <- 10000

ecm_aleat <- numeric(length(n_vals))
ecm_clas <- numeric(length(n_vals))

for (i in seq_along(n_vals)) {
  n <- n_vals[i]
  
  estim_aleat <- numeric(Nrep)
  estim_clas <- numeric(Nrep)
  
  for (iter in 1:Nrep) {
    
    muestra_aleat <- rbinom(n, 1, p_y_aleat)
    media_aleat <- mean(muestra_aleat)
    estim_aleat[iter] <- min(1, max(0, 1.5 * (media_aleat - 1/6)))
    
    
    muestra_clas <- rbinom(n, 1, p_y_clas)
    estim_clas[iter] <- mean(muestra_clas)
  }
  
  ecm_aleat[i] <- mean((estim_aleat - tita)^2)
  ecm_clas[i] <- mean((estim_clas - tita)^2)
}


tabla_resultados <- data.frame(
  n = n_vals,
  ECM_Clasico = ecm_clas,
  ECM_Aleat = ecm_aleat
)

print(tabla_resultados) 
##      n ECM_Clasico    ECM_Aleat
## 1   10 0.020782000 0.0415575000
## 2  100 0.008346730 0.0050118025
## 3 1000 0.007087487 0.0005029204
## 4 2000 0.006994958 0.0002563438
## 5 3000 0.007003293 0.0001653930

2h

B <- 100000
y <- c(rep(1, 20), rep(0, 80))
alpha <- 0.05
z_alpha_2 <- qnorm(1 - alpha/2)

theta_hat <- function(sample) {
  p_hat <- mean(sample)
  theta <- min(1, max(0, 3/2 * (p_hat - 1/6)))
  return(theta)
}

theta_bootstrap <- replicate(B, {
  sample_y <- sample(y, size = 100, replace = TRUE)
  theta_hat(sample_y)
})


IC_bootstrap <- quantile(theta_bootstrap, probs = c(alpha/2, 1 - alpha/2))

sd_boot <- sd(theta_bootstrap)
cota_inf <- theta_hat(y) - z_alpha_2 * sd_boot
cota_sup <- theta_hat(y) + z_alpha_2 * sd_boot
IC_bootstrap_2 <- c(max(0,cota_inf), min(1,cota_sup))
##  2.5% 97.5% 
##  0.00  0.17
## [1] 0.00000 0.14908

2i

calcular_nivel_empirico <- function(theta0 = 1/4, 
                                    alpha = 0.05, 
                                    B = 10000, 
                                    n_vals = c(200, 300, 400, 500, 1000, 1500, 2000)) {
  
  crit_val <- qchisq(1 - alpha, df = 1)
  p0 <- 1/6 + (2/3) * theta0
  
  computar_T <- function(y) {
    n <- length(y)
    S <- sum(y)
    ybar <- S / n
    p_hat <- min(5/6, max(1/6, ybar))
    T <- - 2 * (S * log(1/3) + (n - S) * log(2/3) - S * log(p_hat) - (n - S) * (log(1 - p_hat)))
    return(T)
  }
  
  niveles_empiricos <- numeric(length(n_vals))
  
  for (i in seq_along(n_vals)) {
    n <- n_vals[i]
    
    # B muestras de tamaño n bajo H0
    muestras <- matrix(rbinom(n * B, size = 1, prob = p0), nrow = B, ncol = n)
    T_vals <- apply(muestras, 1, computar_T)
    
    niveles_empiricos[i] <- mean(T_vals > crit_val)
  }
  
  return(data.frame(n = n_vals, nivel_empirico = niveles_empiricos))
}

# n chicos
n_pequenos <- 1:150
resultado_pequenos <- calcular_nivel_empirico(n_vals = n_pequenos)

# n grandes
n_grandes <- c(200, 300, 400, 500, 600, 700, 800, 900, 1000, 1500, 2000, 2500, 3000)
resultado_grandes <- calcular_nivel_empirico(n_vals = n_grandes)

n_comb <- c(n_pequenos, n_grandes)
nivel_comb <- calcular_nivel_empirico(n_vals = n_comb)

Latex

Interactivo

En estos gráfico se observa que:

Este resultado muestra que, si bien el test tiene un comportamiento deficiente para muestras pequeñas, su nivel de significación se ajusta al valor teórico a medida que aumenta el tamaño muestral. Por lo tanto, es recomendable utilizar este test con \(n\) suficientemente grande.

2j

n <- 100
B <- 10000
alpha <- 0.05
theta0 <- 1/4
p0 <- 1/6 + (2/3) * (theta0)
c_crit <- qchisq(1 - alpha, df = 1)

p_vals <- sort(c(seq(1/6, 5/6, by = 0.01), 1/3))
potencia <- numeric(length(p_vals))

calc_aux_logver <- function(p, ybar) {
  ybar * log(p) + (1 - ybar) * log(1 - p)
}


for (i in seq_along(p_vals)) {
  p <- p_vals[i]
  
  stats <- replicate(B, {
    y <- rbinom(n, 1, p)
    ybar <- mean(y)
    
    p_hat <- min(5/6, max(1/6, ybar))
    
    stat <- - 2 * n * (calc_aux_logver(p0, ybar) - calc_aux_logver(p_hat, ybar))
    stat
  })
  
  potencia[i] <- mean(stats > c_crit)
}

Latex

Interactivo